MODULE VEXCSS_MOD
CONTAINS
SUBROUTINE VEXCSS(KIDIA,KFDIA,KLON,PTMST,PRVDIFTS,&
 & PUMLEV,PVMLEV  ,PTMLEV   ,PQMLEV, &
 & PAPHMS,PGEOMLEV,PCPTGZLEV,PCPTS , &
 & PQSAM ,PZ0MM   ,PZ0HM    ,PZ0QM ,PBUOM, &
 & PUCURR,PVCURR  , &
 & YDCST ,YDEXC   , &
 & PCFM  ,PCFH    ,PCFQ )  

USE PARKIND1  , ONLY : JPIM, JPRB
USE YOMHOOK   , ONLY : LHOOK, DR_HOOK, JPHOOK
USE YOS_THF   , ONLY : RVTMP2
USE YOS_EXCS  , ONLY : RLPDD, RLPCC, RLPBB
USE YOS_CST   , ONLY : TCST
USE YOS_EXC   , ONLY : TEXC

! (C) Copyright 1990- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

!     ------------------------------------------------------------------

!**   *VEXCSS* - DETERMINES THE EXCHANGE COEFFICIENTS BETWEEN THE
!                 SURFACE AND THE LOWEST MODEL LEVEL WITH HELP OF
!                 STABILITY AS FUNCTION OF OBUKHOV-L.

!     Original A.C.M. BELJAARS       E.C.M.W.F.    26/03/90.
!     Modified A.C.M. BELJAARS   26/03/99 Tiling of the land surface.
!     Modified J. HAGUE          13/01/03 MASS Vector Functions 
!     Duplicated P. LOPEZ        03/06/05 New version for linearized physics
!     Modified   P. Viterbo ECMWF 12/05/2005 Externalize SURF
!                                   (based on vdfexcss)
!                M. Janiskova    27/06/2005 removed option for MASS vector
!                                   functions not use in corresponding TL/AD
!                M. Janiskova    15/02/2006 code re-organization for efficient
!                                   computation of its TL/AD versions
!                P. Lopez        July 2025  Added ocean currents

!     PURPOSE
!     -------

!     DETERMINE EXCHANGE COEFFICIENTS BETWEEN THE SURFACE AND THE
!     LOWEST MODEL LEVEL

!     INTERFACE
!     ---------

!     *VEXCSS* IS CALLED BY *SURFEXCDRIVERS*

!     INPUT PARAMETERS (INTEGER):

!     *KIDIA*        START POINT
!     *KFDIA*        END POINT
!     *KLON*         NUMBER OF GRID POINTS PER PACKET

!     INPUT PARAMETERS (REAL):

!     *PTMST*        TIME STEP
!     *PRVDIFTS*     Semi-implicit factor for vertical diffusion discretization
!     *PUMLEV*       X-VELOCITY COMPONENT AT T-1, lowest model level
!     *PVMLEV*       Y-VELOCITY COMPONENT AT T-1, lowest model level
!     *PTMLEV*       TEMPERATURE AT T-1, lowest model level
!     *PQMLEV*       SPECIFIC HUMUDITY AT T-1, lowest model level
!     *PAPHMS*       PRESSURE AT T-1, surface
!     *PGEOMLEV*     GEOPOTENTIAL AT T-1, lowest model level
!     *PCPTGZLEV*    DRY STATIC ENERGY, lowest model level
!     *PCPTS*        DRY STATIC ENERGY AT THE SURFACE
!     *PQSAM*        SPECIFIC HUMIDITY AT THE SURFACE
!     *PZ0MM*        AERODYNAMIC ROUGHNESS LENGTH
!     *PZ0HM*        ROUGHNESS LENGTH FOR TEMPERATURE
!     *PZ0QM*        ROUGHNESS LENGTH FOR MOISTURE
!     *PBUOM*        BUOYANCE FLUX AT THE SURFACE
!     *PUCURR*       OCEAN CURRENT U-COMPONENT
!     *PVCURR*       OCEAN CURRENT V-COMPONENT

!     OUTPUT PARAMETERS (REAL):

!     *PCFM*         PROP. TO EXCH. COEFF. FOR MOMENTUM (C-STAR IN DOC.)
!                    AT LOWEST MODEL LEVEL; CALLED WITH ZCFM(1,KLEV)
!     *PCFH*         PROP. TO EXCH. COEFF. FOR HEAT     (C-STAR IN DOC.)
!                    AT LOWEST MODEL LEVEL; CALLED WITH ZCFM(1,KLEV)
!     *PCFQ*         PROP. TO EXCH. COEFF. FOR MOISTURE (C-STAR IN DOC.)
!                    (THIS ARRAY APPLIES TO BOTTOM LAYER ONLY)

!     METHOD
!     ------

!     THE ALGEBRAIC RELATION BETWEEN Z/L AND THE RICHARDSON NUMBER
!     IS SOLVED ITERATIVELY. THE STABILITY FUNCTIONS ARE THE SO-CALLED
!     PROFILE PSI-FUNCTIONS.
!     THE INITIAL GUESS (E.G. FROM PREVIOUS TIMESTEP) IS BACK-
!     SUBSTITUTED TO OBTAIN A SECOND APPROXIMATION. FURTHER ITERATION
!     IS DONE BY LINEAR INTER(EXTRA)POLATION (NEWTON'S  METHOD WITH
!     THE DERIVATIVE FROM SUCCESSIVE APPROXIMATIONS). 

!     ------------------------------------------------------------------

IMPLICIT NONE

INTEGER(KIND=JPIM),INTENT(IN)    :: KLON 
INTEGER(KIND=JPIM),INTENT(IN)    :: KIDIA 
INTEGER(KIND=JPIM),INTENT(IN)    :: KFDIA 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTMST 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PRVDIFTS
REAL(KIND=JPRB)   ,INTENT(IN)    :: PUMLEV(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PVMLEV(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTMLEV(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PQMLEV(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PAPHMS(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PGEOMLEV(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCPTGZLEV(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCPTS(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PQSAM(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PZ0MM(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PZ0HM(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PZ0QM(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PBUOM(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PUCURR(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PVCURR(:) 
TYPE(TCST)        ,INTENT(IN)    :: YDCST
TYPE(TEXC)        ,INTENT(IN)    :: YDEXC
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PCFM(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PCFH(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PCFQ(:)  

!*            LOCAL STORAGE
!             ----- -------
INTEGER(KIND=JPIM) :: JLEN

REAL(KIND=JPRB) ::  Z1DZ0M(KLON)  ,Z1DZ0H(KLON),&
 & Z1DZ0Q(KLON)  ,ZXLNM(KLON),&
 & ZXLNH(KLON)   ,ZXLNQ(KLON),&
 & ZDU2(KLON)    ,ZNLEV(KLON),&
 & ZRICLS(KLON), ZQSAM(KLON) 

INTEGER(KIND=JPIM) :: JIT, JL

REAL(KIND=JPRB) :: Z2B, Z3B, ZAUX1, ZAUX2, ZC, ZCDNH,&
 & ZCDNM, ZCDNQ, ZCH, ZCM, ZCON1, ZCON2, ZCONS12, &
 & ZCQ, ZD, ZDRORO, ZIPBL, ZTPFAC1, ZUABS, ZWST2, ZRGI, ZRKAP3B
REAL(KIND=JPRB) :: ZRIMIN, Z2S, Z5S, Z6S, Z7S
REAL(KIND=JPRB) :: ZDIV1, ZDIV2, ZDIV3, ZDIV4, ZDIV5, ZDIV6
REAL(KIND=JPRB) :: ZDIV7, ZDIV8, ZDIV9, ZDIV10, ZDIV11
 
REAL(KIND=JPHOOK) :: ZHOOK_HANDLE

!     ------------------------------------------------------------------

!*       1.   INITIALIZE CONSTANTS
!             ---------- ----------

IF (LHOOK) CALL DR_HOOK('VEXCSS_MOD:VEXCSS',0,ZHOOK_HANDLE)
ASSOCIATE(RD=>YDCST%RD, RETV=>YDCST%RETV, RG=>YDCST%RG, &
 & REPDU2=>YDEXC%REPDU2, RKAP=>YDEXC%RKAP, RPARZI=>YDEXC%RPARZI)

ZTPFAC1=PRVDIFTS

ZCONS12=ZTPFAC1*PTMST*RG/RD
ZCON1  =RVTMP2-RETV
ZCON2  =2.0_JPRB/3._JPRB
ZRIMIN = 1.E-3_JPRB
ZRGI = 1.0_JPRB/RG

!     CONSTANTS FOR THE STABILITY FUNCTIONS

Z2B=2.0_JPRB*RLPBB
ZD=RLPDD
ZC=RLPCC
Z3B=3._JPRB*RLPBB
ZRKAP3B = RKAP*RKAP*Z3B

!     PBL HEIGHT FOR W* - EFFECT

ZIPBL=RPARZI

!     ------------------------------------------------------------------

!*       2.   COMPUTATION OF BASIC QUANTITIES
!             ----------- -- ----- ----------

DO JL=KIDIA,KFDIA

!             (W*)**2, WIND SHEAR, RICHARDSON NUMBER,

  IF (PBUOM(JL)  <=  0.0_JPRB) THEN
    ZWST2=0.0_JPRB
  ELSE
    ZWST2=(PBUOM(JL)*ZIPBL)**ZCON2
  ENDIF
  ZDU2(JL) = MAX(REPDU2,(PUMLEV(JL)-PUCURR(JL))**2 + (PVMLEV(JL)-PVCURR(JL))**2 + ZWST2)
  ZDIV1 = 1.0_JPRB/(PCPTGZLEV(JL)+PCPTS(JL)-PGEOMLEV(JL))
  IF (PQSAM(JL) < 0.1_JPRB) THEN
    ZQSAM(JL) = PQSAM (JL)
  ELSE
    ZQSAM(JL) = 0.1_JPRB
  ENDIF
  ZDRORO=2.0_JPRB*(PCPTGZLEV(JL)-PCPTS(JL))*ZDIV1 &
   & - ZCON1*(PQMLEV(JL)-ZQSAM(JL))
  ZDIV2 = 1.0_JPRB/ZDU2(JL)
  ZRICLS(JL)=PGEOMLEV(JL)*ZDRORO*ZDIV2

!             COMMON FACTORS IN NEUTRAL FORMULAE AND
!             DRAG COEFFICIENTS.

  ZNLEV (JL)=PGEOMLEV(JL)*ZRGI+PZ0MM(JL)
  ZDIV3 = 1.0_JPRB/PZ0MM(JL)
  Z1DZ0M(JL)=ZNLEV(JL)*ZDIV3
  ZDIV4 = 1.0_JPRB/PZ0HM(JL)
  Z1DZ0H(JL)=ZNLEV(JL)*ZDIV4
  ZDIV5 = 1.0_JPRB/PZ0QM(JL)
  Z1DZ0Q(JL)=ZNLEV(JL)*ZDIV5
  ZXLNM(JL) =LOG(Z1DZ0M(JL))
  ZXLNH(JL) =LOG(Z1DZ0H(JL))
  ZXLNQ(JL) =LOG(Z1DZ0Q(JL))

!*          EXCHANGE COEFFIENTS DEPENDING UPON RICHARDSON NUMBER

  ZCDNH =1.0_JPRB/(ZXLNM(JL)*ZXLNH(JL))
  ZCDNQ =1.0_JPRB/(ZXLNM(JL)*ZXLNQ(JL))
  ZCDNM =1.0_JPRB/(ZXLNM(JL)*ZXLNM(JL))

  ZCM = ZC*SQRT(Z1DZ0M(JL))
  ZCH = ZC*SQRT(Z1DZ0H(JL))
  ZCQ = ZC*SQRT(Z1DZ0Q(JL))

!*        STABILITY FUNCTIONS 

  IF (ZRICLS(JL) > 0) THEN
    Z2S = SQRT(1.0_JPRB+ZD*ZRICLS(JL))
    ZDIV6 = 1.0_JPRB/(Z2S+Z2B*ZRICLS(JL))
    PCFM(JL) = ZCDNM*Z2S*ZDIV6
    ZDIV7 = 1.0_JPRB/(1.0_JPRB+Z3B*ZRICLS(JL)*Z2S)
    PCFH(JL) = ZCDNH*ZDIV7
    PCFQ(JL) = ZCDNQ*ZDIV7
  ELSE
    IF ( -ZRICLS(JL) > ZRIMIN) THEN
      Z2S = SQRT(-ZRICLS(JL))
    ELSE
      Z2S = SQRT(ZRIMIN)
    ENDIF
    Z5S = ZRKAP3B*ZCDNM*ZCM
    ZDIV8 = 1.0_JPRB/(1.0_JPRB+Z5S*Z2S) 
    PCFM(JL) = ZCDNM*(1.0_JPRB-(Z2B*ZRICLS(JL))*ZDIV8)
    Z6S = ZRKAP3B*ZCDNH*ZCH
    ZDIV9 = 1.0_JPRB/(1.0_JPRB+Z6S*Z2S)
    PCFH(JL) = ZCDNH*(1.0_JPRB-(Z3B*ZRICLS(JL))*ZDIV9)
    Z7S = ZRKAP3B*ZCDNQ*ZCQ
    ZDIV10 = 1.0_JPRB/(1.0_JPRB+Z7S*Z2S)
    PCFQ(JL) = ZCDNQ*(1.0_JPRB-(Z3B*ZRICLS(JL))*ZDIV10)
  ENDIF

  ZDIV11 = 1.0_JPRB/(PTMLEV(JL)*(1.0_JPRB+RETV*PQMLEV(JL)))
  ZAUX1  = (ZCONS12*PAPHMS(JL))*ZDIV11
  ZUABS  =SQRT(ZDU2(JL))
  ZAUX2  =ZAUX1*ZUABS*RKAP**2

  PCFM(JL)=ZAUX2*PCFM(JL)
  PCFH(JL)=ZAUX2*PCFH(JL)
  PCFQ(JL)=ZAUX2*PCFQ(JL)

ENDDO

END ASSOCIATE
IF (LHOOK) CALL DR_HOOK('VEXCSS_MOD:VEXCSS',1,ZHOOK_HANDLE)
END SUBROUTINE VEXCSS
END MODULE VEXCSS_MOD
